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Abstract 

Using random matrix technique we determine an exact relation between the eigen- 
value spectrum of the covariance matrix and of its estimator. This relation can be 
used in practice to compute eigenvalue invariants of the covariance (correlation) 
matrix. Results can be applied in various problems where one experimentally esti- 
mates correlations in a system with many degrees of freedom, like for instance those 
in statistical physics, lattice measurements of field theory, genetics, quantitative 
finance and other applications of multivariate statistics. 
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Statistical systems with many degrees of freedom appear in numerous re- 
search areas. One of the most fundamental issues in studies of such systems is 
the determination of correlations. In practice, one encounters frequently the 
following situation: one samples the system many times by carrying out inde- 
pendent measurements. For each sample one estimates values of the elements 
of the covariance matrix, and then takes the average over a set of samples. 
The statistical uncertainty of the average of individual elements of the ma- 
trix generically decreases with the number of independent measurements T 
as ~ 1/yT. There are N(N + l)/2 independent elements of the correlation 
matrix for a system with N degrees of freedom. Thus, naively, the total uncer- 
tainty encoded in the correlation matrix may be expected to be proportional 
to iV(iV + l)/2 and 1/y/T, and therefore to be large for 'non-local' quantities 
which depend on many elements of the correlation matrix even if one performs 
a large number of measurements, of the order of the number of degrees of free- 
dom in the system. It turns out that this naive expectation is far from true. 
Such 'non-local' quantities occur, in particular, in the eigenvalue analysis of 
the correlation matrix [1]. A question which we address in this paper is how 
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the spectrum of the experimentally measured covariance (two-point correla- 
tion) matrix is related to the spectrum of the genuine correlation matrix for 
the system. 



To be specific, consider a statistical system consisting of iV real degrees of 
freedom x iy i = 1, . . . , N with a stationary probability distribution: 



N 



p(xi, . . .,x N ) Yl dx n (1) 



n=l 

such that 

N 



/ Xi p(xi, . . . , x N ) Yl dx n = Mi . (2) 

J n=l 

The covariance matrix for the system is defined as 

N 

Qj = / XiXj p(xi, . . . x N ) Y[ dx n . (3) 

J n=l 

Further, assume that the system belongs to the Gaussian universality class. 
Under this assumption the probability distribution can be approximated by 

N 1-V2 / 1 \ N 

p(x 1 , . . . , x N ) Y[ dx n = [(27r) 7V detCj exp -- J2 XiC^xj J| dx n (4) 

n=l y ij J n=l 

where Cjj is a covariance matrix (3) of the system. By construction it is a 
symmetric, positive-definite matrix. In fact, for a wide class of models, the 
Gaussian approximation well describes the large N behavior of the system 
as a consequence of the central limit theorem. Deviations from the Gaussian 
behavior can result either from the presence of fat (heavy) tails in the prob- 
ability distribution or from collective excitations of many degrees of freedom. 
None of these effects will be discussed here. 



Experimentally, the correlation matrix is computed as follows. One performs a 
series of T independent measurements. Assume T > N. The measured values 
x n form a rectangular N x T matrix X with elements X nt , where X nt is the 
measured value of the n th degree of freedom x n in the t th experiment t = 
1, . . . , T. The experimental correlation matrix is computed using the following 
estimator 

%4EM't = ^ XX % (5) 
1 t=i 1 

where X T is the transpose of X. We expect that for T — > oo the estimated 
values Cij will approach the elements C^. More precisely, if the measurements 
are independent, the probability distribution of measuring a matrix X of values 
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X nt is a product of probabilities for individual measurements 



T / N \ 

p(x)dx = n [p(Xit, ■ ..x m ) n dx nt 

t=l \ n=l / 



(6) 



where 

N,T 

DX = J] d*nt • (7) 

n,t=l 

In particular, for the Gaussian approximation 



/ 1 T \ N,T 

P(X)DX=Afe W [--J2XitC- j 1 X jt ) [] dX nt 

\ 1 t=l I n,t=l 

= ATexp (~7p X T C _1 X^ DX 



where M is a normalization factor which ensures that / P(X)DX = 1. In this 
particular case J\f = [(27r) Ar detC]~ T//2 . All averages over measured values X nt 
are calculated with this probability measure. We shall denote these averages 
by (...). In particular we see that 

(X it Xj t i) = Cij5 u > ■ (9) 

This relation reflects the assumed absence of correlations between measure- 
ments. In general, if measurements are correlated, the right-hand side of the 
last equation can be expressed by a matrix C it jt' in double indices. 

After these introductory remarks, we come back to the problem of relating the 
spectrum of the covariance matrix C to the spectrum of its estimator c. We 
denote the eigenvalues of the matrix C by A n n = 1, . . . , N. For a given set of 
eigenvalues we can calculate matrix invariants, like for example the spectral 
moments 

M k = llV C k = 1 £ A£ = JdA Po (A)A k (10) 
where the density of eigenvalues po(A) is defined as 

Po(A) = ^E5(A-A„) . (11) 

JV 71=1 

The question is how these quantities are related to the analogous quantities 
defined for the estimator of the correlation matrix c 

m fc = l(Trc fe )= jd\p(\)\ k (12) 
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where the eigenvalue density of the matrix estimator is 

pw4(e/( a - a «)} • ( l3 ) 

We expect that the dependence of the estimated spectrum p(A) and the gen- 
uine spectrum po(A) should be controlled by T and N. Indeed, as we shall 
see later, it turns out that for iV — > oo this dependence is governed by the 
parameter r = N/T, which we assume to be finite. 

In order to derive the relation between the spectral properties of the covariance 
matrix and its estimator it is convenient to define resolvents: 

G(Z) = (Z1 N -Cy 1 (14) 

and 

g(z) = ((zl N - c)~^ = ((zl N - ^XXf 1 ^ (15) 

where Z and z are complex variables. The symbol Ijv stands for the N x N 
unit matrix. Expanding the resolvents in \ jZ (or 1/z) one sees that they can 
be interpreted as generating functions for the moments 

1 00 1 

M(Z) = -Tt[ZG(Z)] - 1 = £ ~ k M k (16) 

k=l 

and 

1 00 1 

m{z) = -Tr[zg(z)} - 1 = £ -m fc . (17) 

iV k=i z 

From the relation between M(Z) and m(z) one can determine the correspond- 
ing relation between the eigenvalue spectra po(A) and p(A). Indeed, taking 
the imaginary part of Tr g(z)/N (and Tr G(Z)/N) for z = A + i0 + (or 
Z = A + i0 + ), where A is real, we can directly calculate the eigenvalue densi- 
ties p(A) (and Po(A)): 

p ( A ) = -Ilm-|-Trg(A + 20 + ) (18) 

71 iV 

as follows from the standard relation for distributions: (x + iO + )~ 1 = PVx^ 1 — 
inS(x), where PV stands for principal value. 

The fundamental relation between the generating functions (16) and (17) is de- 
rived in the Appendix by means of a diagrammatic technique [3] for calculating 
integrals (15) with the Gaussian measure (8). The large N limit corresponds to 
the planar limit in which only planar diagrams contribute. This significantly 
simplifies considerations and allows one to write down closed formulae for the 
resolvents. 
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This fundamental relation between the generating functions (16) and (17) 
reads 

m(z) = M{Z) (19) 
where the complex number Z is related to z by the conformal map (56): 

(20) 



1 + rm{z) 

or equivalently, if we invert the last relation for z = z(Z), as: 

z = Z(l + rM(Z)). (21) 

The equations (19,20) were already announced in our earlier work [2]. They can 
for example be used to compute moments of the genuine correlation function 
C from the experimentally measured moments of the estimator c. Indeed, 
combining (19) and (20) we obtain the following equation: 

m(z) = M (- Z —- j (22) 

v ' \l + rm(z) J v ' 

which gives a compact relation between moments and M k . 



k=l Z k=l Z \ 1=1 Z / 

from which we can recursively express by Mi, I — 

m 1 = Mi 

m 2 = M 2 + rMl 

m 3 = M 3 + :5r. U ; .U, + r 2 M\ 

or inversely: M k by m h I = 1, . . . , k: 

Mi = mi 

M 2 = 1712 — rm\ 

M 3 = m 3 — 3rmim 2 + 2r 2 m\ 



(24) 



(25) 



Let us observe that for r < 1 the functions M(Z) and m(z) can also be 
expanded around z = Z = 0. In this case 

oo 

M(Z) =-J2 Z h M_ k , (26) 

k=0 
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where 

M_ fc = ^TrC~ k . (27) 

Similarly 

oo 

™( z ) = -J2 zkm -k, (28) 

fc=0 

where 

m _ k = ±(Trc- k ). (29) 
Using the same manipulation as before we obtain 

oo oo oo 

]T .1/ ,Z ; = ]T m_ fc Z fe (l -r-rX; M_,Z')* (30) 
fc=i fc=i j=i 

and hence: 

M_! = (1 - r)m_i 

M_ 2 = (1 — r) 2 m_ 2 — r(l — r)m 2 , 

V ) \ ) -i (31) 

M_ 3 = (1 — r) 3 m_ 3 — r(l — r) 2 m_ 1 m_ 2 — r 2 (l — r)m^ x 



The relations between moments can be used directly in practical applications 
to clean the spectrum of the correlation matrix. Before discussing this let us 
make a comment. The formulae (19) and (20) encode full information about 
the relation between the eigenvalue spectrum po(A) and p(A) for a given r. 
In particular, if one knows the spectrum po(A) of the correlation matrix C 
one can exactly determine for a given r the shape of the spectrum p(A) of the 
estimator dressed by statistical fluctuations. One does it as follows. From the 
eigenvalue spectrum po(A) one deduces an explicit form of the function M(Z) 
and of the right hand side of the equation (21). Inverting the equation (21) for 
Z one finds the dependence Z = Z(z). Inserting it to the equation (19) one 
determines the function m(z). Taking the imaginary part along the cuts of the 
map m(z) on the real axis (18) one eventually finds p(X). One can easily write 
a numerical program which realizes this procedure. In few cases the solution 
is possible analytically. Let us shortly discuss them. 

Consider the correlation matrix C whose spectrum is given by a sequence of 
degenerate eigenvalues p,i, i — 1, . . . , K with degeneracies rii. Consequently, 
defining pi = rii/N, J2iPi — 1; we have 

M(Z) = f; f^- . (32) 

i=l & — 

This form is particularly simple to discuss. One should however keep in mind 
that the relations (19,20) remain valid also in a more general case, for instance, 
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when in the limit N — > oo the spectrum of po(A) is not a sum of delta functions 
but approaches some continuous distribution. The map (21) now reads 

(33) 

Clearly, if we solve this equation for Z = Z(z) we obtain a multi- valued 
function, except in the case r = for which we have a simple relation z = Z. 
The "physical" Riemann sheet of the map Z = Z(z) is singled out by the 
condition Z — > z for z — > oo. On this sheet the complex z-plane is mapped on a 
part of the Z plane without a simply or multiply connected region surrounding 
the poles at Z — Hi. 

As an illustration, let us consider the simplest case, where K — 1. In this case 
we have only one eigenvalue \i\ = \i and p± — 1 and M(Z) = fJ>/(Z — ji). The 
map (20) has the following form 

z = Z + r^^. (34) 
Zi — n 

If one rewrites the right-hand side of this equation using polar coordinates 
(R, 0) around the pole: Z — /i = Re 1 ^: 

2 

z = R e i<t> + nL e -i4> + II (i + r ) (35) 
R 

one can see that the equation is invariant under the "duality" transformation: 

2 

R < — > —— , < — > —(p (36) 
R 

which maps the inside of the circle \Z — fi\ = [i^fr onto the outside and vice 
versa. Obviously, the outside corresponds to the "physical" Riemann sheet of 
the inverse map Z — Z(z) 

Z= l -({l-r)p + z + ^/(z -//+)(* -//_)) , (37) 

since in this region Z ~ z for z — > oo. The two constants in the last equation 
are [x± = ± \fr) 2 . Along the cut /x_ < z < fi + on the real axis the map 
Z = Z(z) becomes complex and ambiguous: it has a phase (sign) ambiguity 
which is related to the fact that the cut is mapped into the limiting circle 
where the two Riemann sheets meet. 

From (19) we can easily find the generating function m(z) and then from (18) 
the spectral density of the correlation matrix c: 

1 \/(u+ - A) (A - a J) 
p(A) = -— ^— ^ — . (38) 



7 



This is a well known result in random matrix theory [4] for the spectral dis- 
tribution of the Wishart ensemble. It is interesting to interpret this result as 
a statistical smearing of the initial spectral density po(A) given by the delta 
function localized at p into a wide peak p(A) supported by the cut 
due to a finite series of measurements. The larger r the larger is the width of 
the resulting distribution p(A). For r = 1 the lower limit of the distribution is 
at /!_ = signaling the appearance of zero modes in the matrix c. One can 
show by considering anti- Wishart matrices that above the limiting value r = 1 
the zero mode sector of c grows with r reflecting an increasing indeterminacy 
of the spectrum of the covariance matrix C when the underlying statistical 
sample becomes too short. 

As a second example we consider the case, where the genuine covariance matrix 
C has two different eigenvalues p±, p 2 with relative weights p±, P2, pi +P2 = 1- 
In this case we can also find an explicit form of the map Z(z) solving the 
corresponding cubic (Cardano) equation. Depending on the parameters pi, pi 
the map Z(z) has one or two cuts on the real axis on the z plane, which means 
that the corresponding eigenvalue distribution p(X) has a support on one or 
two intervals (fig.l). It is a simple exercise to find the critical value r c of r at 
which a single cut solution splits into a two-cut one: 

- 3. (39) 



' ((Pi/i?) 1/3 + MY^Y 

For instance, this formula gives r c as 0.01 if the correlation matrix C has two 
eigenvalues p± = 1 and p 2 = 1.1 and the corresponding weights p\ — p 2 — 1/2. 
Thus in this case, to observe a bimodal signal in the measured spectrum one 
has to perform T measurements with T of order 100 N. If the gap Ap between 
the eigenvalues p 1 and p 2 is larger the binomial structure is visible already for 
smaller T. As an example we show in fig.l typical shapes of the distribution 
for pi — 1, p 2 — 2, pi = p 2 = 1/2 for r < r c , r = r c and r > r c . The shape 
of the spectral functions for r in the range from 1 to, say, 0.2 — 0.3 resembles 
that for a single eigenvalue. When r approaches r c it starts to deviate from 
this shape developing a double peak structure which eventually splits into two 
separate parts at r = r c . If r is further decreased the two peaks get narrower. 
They eventually entirely localize at p\ — 1 and p 2 = 2 for r = 0, that is for 
an infinitely long sample for which the spectrum of the underlying correlation 
matrix is recovered. 



The method can be straightforwardly generalized from K — 1,2 to arbitrary 
K, pi, . . . , px with J^Pi — 1, although only the K = 3 case is solvable analyt- 
ically (quartic Ferrari equation). In other cases one can use a numerical imple- 
mentation of the general procedure, which we described before, to determine 
the shape of the spectrum of the estimator p(X) from any given distribution 
Po(A) and for any r. 



8 



r = 0.01 




0.5 1 1.5 2 2.5 3 3.5 4 4.5 

A 

Fig. 1. The figures represent spectra of the eigenvalue distributions p(A) of 
the experimental correlation matrix measured in a series of measurements for 
r = 0.01,0.115,0.3, respectively. The underlying correlation matrix has two eigen- 
values jUi = 1 and fi2 = 2 with the weights Pi = P2 = 0.5. At the critical value 
r = r c = 0.115 (eq. 39) the spectrum splits. The spectral densities are calculated 
analytically. 

In practice one is however interested in the opposite problem that is in the 
determination of the spectrum po(A) of the genuine correlation matrix C from 
the distribution of the measured eigenvalues. This is a difficult problem for 
the following reason. Having one sample of T measurements of a system with 
./V degrees of freedom one has as a result only N eigenvalues of the correlation 
matrix. From N random numbers it is impossible to reconstruct accurately 
the exact form of the underlying distribution p(A) according to which they 
are distributed. Thus the input function p(A) for the procedure leading from 
p(A) to po(A) has a large statistical uncertainty. In effect the output function 
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Po(A) is not well controlled. 



In some applications one is not interested in the exact shape of the spectral 
function po(A) but in its moments. In such cases one can directly make use 
of the relations (25) between the measured moments m& and the moments of 
the original distribution po(A) to determine the latter ones. 

Moreover, the relations (25) between the moments may in some cases signif- 
icantly improve the determination of the eigenvalues A of the original dis- 
tribution. Assume that on top of the purely statistical information obtained 
by independent measurements of a system consisting of iV-degrees of freedom 
we have at our disposal some additional non-statistical knowledge about the 
system. For example, assume that we know that degrees of freedom can be 
grouped into K sectors of degenerate independent constituents. Each sector 
is represented by an eigenvalue pk and a degeneracy rik, or equivalently by 
the fraction pk = rik/N of eigenvectors related to this eigenvalue. In other 
words, we assume that the resolvent M(Z) (16) of the correlation function C 
is given by the formula (32) with a specific K. Thus the problem of determining 
eigenvalues Aj is reduced to the problem of determining parameters pk, Pk- If 
K <C N the problem has much less unknowns. This non-statistical knowledge 
of the system can be used as follows. We can explicitly express the moments 
Mfc by the yet unknown parameters pj, pj. We denote the corresponding func- 
tions by Ml h (pj, pj). On the other hand we can measure experimentally the 
moments and from the relations (25) we can find the corresponding values 
which we denote by Ml xp {rrij). Using the jack-knife procedure we can also 
estimate the statistical errors of M% xp, s. Minimizing x 2 (Pj,Pj) : 



we can eventually find optimal values of the parameters Pj,Pj- Let us make 
a few comments. Obviously, L must be equal or larger than the number of 
unknown parameters. If it is equal then the minimization of x 2 amounts to 
solving equations Ml h (pj,pj) = M^ xp . In practice, if possible, L should be 
taken larger than the number of free parameters. In this case the weights 
1/A| in x 2 take care of the gradually decreasing importance of the higher 
moments which are usually estimated with larger errors. 

Let us illustrate the method at work using the following exercise as an example. 
The exercise has two parts. First we generate T x N matrix X from the 
Gaussian distribution (8) with a given covariance matrix C which has K 
eigenvalues pk with the weights Pk, k — 1, . . . , K. Then we treat the parameters 
Pk,Pk as unknown, and experimentally determine their optimal values using 
the method of moments. In the end we compare the measured values with those 
used in the generator. We repeat this procedure many times to estimate the 




(40) 
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0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 

A 

3 | 1 1 1 1 




0.5 1 1.5 2 2.5 

A 



Fig. 2. The figures represent eigenvalue distributions measured in n = 10 5 measure- 
ments for N = 100, T = 333 (upper); the cleaned spectrum using the additional 
information that there are exactly two eigenvalues (middle); and using the infor- 
mation that there are two equally probable eigenvalues in the correlation matrix 
(lower). The peaks in the last two figures localize around correct values fJ-i = 1 and 
/i2 = 2 despite the number of measurements is relatively small: N/T = r = 0.3. 

error one makes, when one estimates the spectrum in a single measurement. 

As a first example we consider the case K — 2, /ii — 1, /i 2 — 2 and p\ = p 2 = 
0.5 as before and with = 100 and r = N/T = 0.3. The theoretical curve of 
the distribution p(\) is shown again in fig. 2 where it is compared with the 
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1 2 3 4 5 6 

A, A 

Fig. 3. The figure represents eigenvalue distributions measured in n = 10 5 measure- 
ments for N = 100, T = 333 (solid line) and the cleaned spectrum using additional 
information that there are three equally probable eigenvalues (dashed line); The 
peaks in the right figure localize around correct values fi± = 1, p2 = 2 and ^3 = 3 
despite the asymmetry parameter is r = 0.3. 

experimentally determined distribution. The experimental curve is obtained 
by averaging over n = 10 5 independent experiments. In each experiment the 
matrix X was generated and the eigenvalues of c were calculated. Thus the 
resulting experimental histogram was constructed out of n ■ N = 10 7 eigenval- 
ues. One sees that the experimental curve fits very well to the theoretically 
predicted. Actually, one cannot distinguish by bare eye the two curves on the 
first plot in fig. 2 . If we knew only this curve we would not be able to con- 
clude from it that the underlying covariance matrix has only two eigenvalues. 
Now assume that we know that there are exactly two eigenvalues and use the 
method of moments to compute them. We have to measure at least L = 3 
moments mi, I — 1, 2, 3 in order to determine by minimizing \ 2 (40) the three 
unknown parameters p\, p 2 , p\. The resulting spectrum po(X) averaged over 
n experiments is shown in the next plot in fig. 2. The spectrum clearly shows 
the double peak structure with the correct localization of peaks around /ii = 1 
and \ii = 2. The areas under the peaks are approximately equal which means 
that pi ~ P2 ~ 0.5, as expected. Assume that on top of the information that 
there are exactly two eigenvalues we additionally know that both are equally 
probable p\ = P2 = 0.5. In effect, we have only two unknown parameters pi 
and fi2- As we see in fig. 2 this additional information leads to the further 
sharpening of the resulting spectrum around the expected values pi — 1 and 
P2 = 2. Another example, for K = 3 and p,\ — 1 p,2 — 2 and ps = 3 is shown 
in fig. 3. The method of moments works very well and allows us to localize 
the positions of eigenvalues of the underlying covariance matrix C. It is clear 
that additional information considerably reduces the error of the estimate. 

To summarize: In the limit N — > oo we give an exact relation between the 
spectrum of the true correlation matrix between the degrees of freedom in a 
statistical system and the averaged spectrum of the natural estimator of this 
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correlation matrix used in a statistical sampling. The latter depends on the 
parameter r = N/T. For finite r the spectrum of the estimator is smeared 
and hides behind the statistical noise the information about the positions of 
eigenvalues of the exact correlation matrix. The averaged spectrum is very 
similar to that of a spectrum obtained from a single matrix Cy (5) used in 
practical applications. The exact relation can be formulated as an infinite 
sequence of linear relations between the moments of the true correlation matrix 
and of its averaged estimator. In principle both sets are equivalent. In practice, 
the estimator is known only as a single realization of the measuring process and 
the estimated moments are known only approximately. The relations between 
the moments can be nevertheless used to determine the spectrum of the true 
correlation matrix from a single estimator. We illustrate the accuracy of the 
method of moments on simple examples. The method of moments can be 
used in practical applications. To give an example, the spectral analysis of 
the covariance matrix plays the central role in the portfolio assessment. The 
knowledge of the covariance matrix is crucial for the optimal asset allocation. 
Typically one constructs the covariance matrix for N of order 500 assets' price 
changes (e.g. SP500 data) and estimates it using a sample of say four years 
of a daily data. In this case T is of order 1000 and hence r = 0.5. As we 
have learned, in this case one almost entirely looses the signal coming from 
the eigenvalues of the genuine covariance matrix. On the other hand, it is 
known that in practice there are only few collective sectors on the market. 
These sectors are represented by K <C iV significant eigenvectors and related 
eigenvalues. One can use the method of moments to localize them [5]. 
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Appendix 

We apply a diagrammatic method [3] to calculate the integral (15). Expanding 
the expression under the average (15) we obtain a series of polynomials in X, 
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which can be integrated term by term for the Gaussian measure ?(X)DX 
(8). It is convenient to write the matrix c as a product of three matrices 
c = X^X r , the first of which is an NxT matrix, the second is TxT , and the 
third TxN. This reveals a sandwich structure of the terms of the expansion: 




The following graphical representation is used in the last equation. The matrix 
X has two types of indices X it , i = 1, . . . , N and t — 1, . . . , T. We shall call 
the first an index of the TV-type, and the second an index of the T-type. 
Indices of the iV-type are drawn as filled circles and of the T-type as empty 
circles, a matrix 1n/z having both indices of the iV-type is drawn as a solid 
line joining two filled circles corresponding to these indices, while a matrix 
1t/T as a dashed line which joins two empty circles representing two indices 
of the T-type of the matrix. An insertion of X is drawn as an ordered pair of 
neighboring circles consisting of a filled circle followed by an open one, while 
an insertion of X r as an ordered pair of an open and a filled circle. 

The X insertions have to be integrated with the Gaussian measure P(X)TJX. 
The terms ^f- and are constant from the point of view of the integration. 
The integration of the X insertions amounts to calculating the correlation 
functions (X iltl . . . X i2kt2k ). Using the Wick theorem we consecutively calculate 
2/c-point correlation functions: the two-point correlation function 

(X it X jt >) = CijS tt ' (42) 

the four-point 

(Aj ltl Xj 2t2 Xj 3t3 Xj 4t4 ) = (X iltl X i2t2 ) (X i3t3 X i4t4 ) 

+ (X iltl X i3t3 ) {X i2t2 X ut4 ) (43) 
+ (X iltl Xi 4t4 ) (Xi 2 t 2 Xi 3 t 3 ) 

and higher correlation functions. All odd correlation functions vanish. All even 
functions can be expressed as sums of all distinct products of two-point func- 
tions given by (42). This observation leads to a particularly simple graphical 
representation: a two-point correlation function is drawn as a double line con- 
sisting of a solid line which connects two filled circles, and of a dashed line 
which connects two empty circles representing indices of X and X r . We asso- 
ciate the contribution with the solid line and 8 t t' with the dashed line. We 
will draw the double lines as semi-circles. Using this convention the four-point 
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correlation function can be drawn as 



(Xi ltl Xi 2t2 Xi 3 t 3 X illt4 ) — 




(44) 



iltl t'2t2 izt'A t^%4 i\t\ t2l2 23*3 tAlA i\t\ t2l2 ^3*3 tA^A 



Putting together equations (41) and (44) we eventually arrive at the diagram- 
matic representation of g(z): 




We are interested in the limit N — > oo, r = N/T = const. In this limit only 
planar graphs contribute to g(z), while non-planar diagrams are suppressed by 
factors proportional to powers of 1 /N — > 0. This property follows from a simple 
counting argument. Each closed internal line of a diagram gives a factor iV (or 
T) coming from taking a trace. Thus a closed solid line contributes a factor 
proportional to N, and a closed dashed line to T. Because we are interested 
in the limit of a constant ratio N/T, T is proportional to N. Each horizontal 
dashed line introduces a factor 1/T ~ 1/N. For a planar diagram the number 
of closed loops is equal to the number of dashed horizontal lines and thus the 
powers of N and of 1/N cancel giving a contribution of order unity, while for a 
non-planar diagram, like for instance the last one displayed in (45), the power 
of 1/N coming from the number of dashed horizontal lines is larger than the 
power of N coming from the number of closed loops, and thus effectively the 
contribution of the diagram has at least one unbalanced factor 1/N which 
makes it vanish in the limit N — > oo, r = N/T = const. 

From here on we neglect the non-planar diagrams. In parallel to the generating 
function g(z) one can define a generating function: 

g *W = ( n T -W ) (46) 

which contains all diagrams of the same shape as in (45) but with dashed 
and solid lines replaced. It is convenient to additionally introduce a class of 
one-line irreducible diagrams having the property that they cannot be divided 
into two separate diagrams by cutting a single horizontal line. Denote the 
sum of one-line irreducible diagrams with two external vertices of the iV-type 
by an d with two external vertices of the T-type by £*(z). The sum 
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of diagrams g(z) can be expressed as a geometric series of the blocks 
connected by solid horizontal lines and similarly the sum of diagrams g*(z) 
by the "E*(z) blocks connected by dashed lines: 



1 Su 1 £u) — S(z) h 

z z z z z z 



1 



/ 1 



JIt-^X/ T1 t -^(z) 
The last equations can be represented graphically as: 



(47) 
(48) 



/ S* \ 

O -D 



■o o- 



-o o- 



-o o--o o - 



The reason why it is useful to introduce the one-line irreducible sums 
and ~E*(z) is that in the large N limit one can write down two additional, 
independent equations which relate £(z) and to g(z) and g*(z). These 

are the following Dyson-Schwinger relations: 



(49) 



which follow from the observation that any planar diagram in the sum g(z) or 
g*(z) can be transformed into a one-line irreducible planar diagram by adding 
an arc joining its external points. 

In this way one obtains a closed set of equations: one has four equations for 
four unknown matrices g, S, g* and £*: 







^ Z) = z^W) 

g * {z) = T-^(z) 
£(z) = CTr[g*(z)' 

£,(*) = l T Tr[g(2)C 



(50) 

(51) 

(52) 
(53) 



One can eliminate X, g* and £* and solve it for g(z) as a function of C. 
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It is convenient to write the result using a variable Z defined by 

zg(z) = ZG(Z) (54) 
where G(Z) is given by equation (14). This amounts to choosing Z as 

z ^)- (55) 

as follows from (50) and (52). The equations (51) and (53) complete the solu- 
tion with the result 

(56) 



l + r(-l + ±Tr[zg(z))) ' 

The equations (54) and (56) relate the resolvent g(z) calculated in a series of 
independent experiments to the genuine resolvent G(z) for the system. 

In the limit of infinitely many measurements T = oo, and iV fixed, which 
corresponds to r = 0, one obtains Z = z and hence g(z) = G(z). As expected, 
in this limit, spectral densities p(X) and po(A) are identical. 
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